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Abstract. Non-linear effects on supernova neutrino oscillations, associated with 
neutrino- neutrino interactions, are known to induce collective flavor transformations 
near the supernova core for ^13 7^ 0. In scenarios with very shallow electron density 
profiles, these transformations have been shown to couple with ordinary matter effects, 
jointly producing spectral distortions both in normal and inverted hierarchy. In this 
work we consider a complementary scenario, characterized by higher electron density, 
as indicated by shock-wave simulations during a few seconds after bounce. In this case, 
early collective flavor transitions are decoupled from later, ordinary matter effects. 
Moreover, such transitions become more amenable to both numerical computations 
and analytical interpretations in inverted hierarchy, while they basically vanish in 
normal hierarchy. We numerically evolve the neutrino density matrix in the region 
relevant for self-interaction effects, using thermal spectra and a representative value 
sin^ 013 = 10^**. In the approximation of averaged intersection angle between neutrino 
trajectories, our simulations neatly show the collective phenomena of synchronization, 
bipolar oscillations, and spectral split, with analytically understandable features, 
as recently discussed in the literature. In the more realistic (but computationally 
demanding) case of non-averaged neutrino trajectories, our simulations do not show 
new significant qualitative features, apart from the smearing of "fine structures" such 
as bipolar nutations. Our results seem to suggest that, at least for non-shallow matter 
density profiles, averaging over neutrino trajectories plays a minor role in the final 
outcome. In this case, the swap of Ve and f^^T- spectra above a critical energy may 
represent an unmistakable signature of the inverted neutrino hierarchy, especially 
for 013 small enough to render further (ordinary or even turbulent) matter effects 
irrelevant. 
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1. Introduction 

Core-collapse supernovae (SN) provide us with an interesting laboratory for studying 
both neutrino properties and their interplay with astrophysical processes (see [1, 2] for 
recent reviews). Mikheev-Smirnov-Wolfenstein (MSW) effects induced on neutrinos by 
background matter [3] have been widely studied as a tool to probe both neutrino masses 
and mixings and the SN dynamics. Recent examples include the characterization, at 
the level of observable SN neutrino event spectra, of shock-wave evolution effects after 
bounce [4-13] and of it possible "erasing" by stochastic matter fluctuations induced by 
turbulence [14-16]. 

Ordinary MSW effects (and their "stochastic" smearing, if any) typically occur 
when a; ~ A, where 

is the vacuum oscillation frequency (in terms of the neutrino energy E and squared mass 
difference Am^)| while 

X{r) = V2GFN,-{r) (2) 

is the Ue-^x interaction energy difference in matter, N^- (r) being the net electron number 
density at the point r. For typical shock- wave density profiles as used, e.g., in [7,8], 
the condition a; ~ A occurs after a few hundred (or even a few thousand) kilometers, 
during the first few seconds after core bounce. Significantly shallower electron density 
profiles (as those adopted in [17] in the context of models with successful r-process 
nucleosynthesis) can instead trigger MSW effects much earlier, around O(IOO) km. 

Besides electrons, neutrinos can also be a nontrivial background to themselves when 
their density is large enough [18, 19]. Self-interaction effects, being inherently non-linear, 
are very different (and much less intuitive) than ordinary MSW effect, and can lead 
to collective flavor transition phenomena in which neutrinos (or antineutrinos) of any 
energy behave similarly [20-34]. The interest of such effects for the neutrino flavor 
evolution in the dense SN core has long been recognized [35-39]. Roughly speaking, 
significant self-interaction phenomena are expected when //(r) >a;, where 

fi{r) = V2GF[N{r) + N{r)] (3) 

N{r) and N{r) being the total effective neutrino (z/g + Ux) and antineutrino {Vg + Vx) 
number density, respectively (to be precisely defined later). 

In the most general case, systems with dense matter and dense neutrino gases are 
thus governed by (at least) three characteristic frequencies: uj (spread over ~ 2 orders 
of magnitude for typical energy spectra); A (roughly decreasing as the third power of 
the distance); and n (decreasing as the fourth power of the distance [17]). Neutrino 
flavor evolution becomes then a complicated multi-scale dynamical problem with a rich 

I In this paper we neglect drn^ = ~ mi <^ AmP, where Am^ = |m§ — 2!- The only relevant 
mixing angle is then ^13, governing the oscillation amplitude in the channels — > and Ve — > Vx 
{x = n or t). 
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phenomenology, which may involve both collective and MSW effects, the latter being 
often assumed (at least in older literature) to lead to flavor transitions. In the SN 
neutrino context, this "old" paradigm (focussing on MSW effects) has dramatically 
changed after the emergence of dominant collective phenomena (such as the so-called 
"bipolar oscillations") studied in detailed, large-scale computer simulations [17,40] as 
well as in simplified but analytical models [41-47] . 

The interaction strength between two (anti)neutrinos is modulated by a factor 
(1 — cos'&ij), where -dij is the angle between their intersecting trajectories. If one 
ignores the spread of , and averages it out along a "representative" radial trajectory 
(single-angle approximation), and if one also assumes that self-interaction effects do not 
interfere with the ordinary MSW ones, then the following picture appears to emerge 
in supernovae from analytical considerations. Nothing relevant occurs for normal 
hierarchy (ms > ^1^2), while, for inverted mass hierarchy (ma < 2), any value of 
^13 7^ (no matter how small [42,43]) can trigger collective pair-conversions of the 
kind UePe t^x^x within the first O(IOO) km. Then, as recently emphasized in [45,46], 
the end of collective effects is marked by a spectral pair-conversion which is complete 
for I7's, while for i/'s it occurs only above a characteristic energy set by lepton number 
conservation. Such spectral "split" (or "stepwise swap" of flavors) would then represent 
an unmistakable signature of sclf-intcraction effects [17, 40, 44-48]. Its robustness needs, 
however, to be better investigated in increasingly reflned SN models including, e.g., 
variable neutrino crossing angles 'dij (multi- angle simulations), which might induce 
kinematical decoherence effects, ("depolarization" and "smearing of oscillations"), which 
are neither obvious nor entirely clear in the few numerical [17, 49] and analytical [34, 43] 
studies performed so far. In our scenario, it turns out that main results are rather robust 
when passing from single- to multi- angle simulations. 

The purpose of this paper is to explore such "self-interaction dominated" scenario 
in a realistic case characterized by: (1) an appropriate matter proflle, where collective 
effects fully develop before MSW effects (if any); (2) continuous, thermal energy 
spectra with significant neutrino-antineutrino (and neutrino flavor) asymmetry, and 
(3) numerical simulations in single- and multi-angle cases. In this sense, our work is 
complementary to the simulations in [17], where the shallower adopted matter profile 
allows MSW effects to occur within (not beyond) the range of collective transitions — 
which leads to a richer phenomenology, but also to more difficult analyses and much 
greater numerical challenges. In the terminology of [42], we study the scenario with 
"thick" rather than "thin" envelope. 

The plan of our work is as follows. In Section 2 we describe our supernova reference 
model. In Section 3 we set the notation and write the neutrino evolution equations 
in single-angle approximation. In Section 4 and 5 we discuss single-angle analytical 
and numerical solutions, respectively. In Section 5 we tackle multi-angle simulations, 
and show that the main final effect (the spectral split) is robust. Conclusions and 
perspectives are presented in Section 6. Technical aspects are discussed in the Appendix. 
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2. Reference Supernova Model 

Our reference model is characterized by the following choices for the initial neutrino 
energy spectra, the geometry and intensity of neutrino emission, and the matter profile, 
which govern the distribution of the three basic "frequencies" cu, /i, and A, respectively. 

We shall use such choices in the usual, simplified context of pure two-family (i/g, Vx) 
evolution [42,43], where Uj. represents a single active flavor. It is worth noticing that, 
for A 7^ 0, this case does not exactly represent the two-family limit (for Sm"^ = 0) of 
the general three-family case, contrary to the familiar MSW effects in SNe. In fact, in 
the presence of self-interactions, the flavor evolution depends on the absolute effective 
neutrino densities and, thus, also on the total number of neutrino families (either two or 
three) assumed to share the total luminosity. The full case (and its proper 2i/ limit) 
will be studied elsewhere. 



2.1. Initial Neutrino Energy Spectra 

We assume normalized thermal spectra with different temperatures T = 1/(3 for Ug, V^, 
and i/x (the latter having the same properties as V^). More precisely, the initial energy 
spectra (j)'^{E) are of the form 

where — 1.202. The average values of E and E'^ are then: 

(E±i) = j dE 4){E) E^^ = c± T±i = c± (5^\ (5) 

where c+ = Ttt'^/ISOCs — 3.151 and c_ = tt^/ISCs — 0.4561. In numerical calculations 
we choose {Ee) = 10 MeV, (Ee) = 15 MeV, and (E^) = (Ex) = 24 MeV for u^, T^e, t^x 
and Vx, respectively, corresponding to 

(3e = 0.315 , pg = 0.210 , (3^^p^^ 0.131 (MeV"^) . (6) 

2.2. Emission Geometry and Intensity 

We adopt the "bulb model" emission described in [17], to which the reader is referred 
for further details. Neutrinos are assumed to be half-isotropically emitted above the 
neutrino-sphere, i.e., all the outward moving modes are equally occupied as expected 
from a blackbody emission. The neutrino number flux per unit energy (in any 
direction) is then given by [17] 

^^^^^ - 2^ W ^ ^ 



where 
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Initial neutrino and antineutrino fluxes 




E (MeV) E (MeV) 

Figure 1. Initial fluxes (at r — 10 km, in arbitrary units) for different neutrino species 
as a function of energy. The fluxes are all proportional to {E) / (E) . 



Rv being the neutrino-sphere radius, while Ly is the total emission power for a given 
neutrino species. In numerical calculations, we assume reference values Ry = 10 km and 
Ly = 10^^ erg/s for each species z/ = u^, V^, I'x, T^x- 

Figure 1 shows the initial neutrino number fluxes per unit energy in arbitrary units 
(all fluxes being proportional to (j)^{E)/{E) through the same normalization constant). 
Notice the significant difference (asymmetry) between neutrinos and antineutrinos, and 
between different neutrino flavors. However, the and fluxes happen to coincide 
at an energy Eeq — 19 MeV, while for the and fluxes the equality occurs at 

~ 24 MeV. Flavor transformations of any kind are not operative for neutrinos at 
E = Eeq, and for antineutrinos a.t E = E^q. 

The spherical symmetry of emission reduces to a cylindrical symmetry along the 
radial line-of-sight (polar axis). At any radius r > Ry along the polar axis, neutrinos will 
arrive with different momenta p characterized by |p| = E, incident polar angle and 
azimuthal angle ip. In the calculation of self-interaction effects, the effective differential 
neutrino number density dup with momentum between p and p + dp is then [17] 

dup = jy{E)dQ = jy{E) dip d cos {} , (9) 

within the cone of sight of the neutrino- sphere, with {} G [0,^9max], being 

^max = arcsin(i?^,/r) . (10) 

In general, angular coordinates are important, since the interaction strength 
between two neutrinos of momenta p and q depends on their relative angle "i^pq through 
the factor (1 — cos?9pq). Calculations embedding the full angular coordinates are dubbed 
"mult i- angle." The often used "single-angle" approximation consists in averaging the 
angular factor along the polar axis, which is assumed to encode the same flavor history 
of any other neutrino direction. In this case, the effective neutrino number density n 
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per unit volume and energy turns out to be [17] 

n{r,E) = 2TiD{r)ju{E) = F^D{r) 



{E) 



(11) 



for the various neutrino species (n = Ug, rig, or n^), where the (species-independent) 
geometrical function D[r) is given by 

2 



D{r) 



1 - 



r 



(12) 



decreasing as for large r. The v-V asymmetry of the model implies rie^ne- 

It is useful to integrate the effective densities per unit energy and volume (rig, rTg, 
and rix — nx) to obtain the effective number densities of i^g, and {yx)-i 



dEn, = F, 



Dir) 



{Ee 



= / dEn 



(13) 
(14) 
(15) 



as well as the total effective number densities of neutrinos and antineutrinos, 

TV = TV, + TV, , (16) 

7V = ]Ve + ]V, , (17) 

from which one can finally derive the parameter ji — \/2Gp (N + N), which governs the 
neutrino self-interaction strength. 

Figure 2 shows the function //(r) in our reference SN model (together with the 
ordinary MSW strength A(r) defined below) in the range r e [10, 200] km relevant for 
self-interaction effects. Also shown are the approximate ranges where these effects induce 
synchronization, bipolar oscillations and spectral split, as discussed later in Sec. 3. 



2.3. Matter Profile 

Shock waves propagating during the first few second after bounce can produce very 
interesting MSW effects when A ~ cj, provided that sin^ ^^13 is not too small [4-11]. 
For the typical (time-dependent) matter profiles used in these studies, such effects 
usually emerge in the first ~ 10 s after a few hundred kilometers from the SN core, in 
which case collective neutrino self-interaction effects are already completely developed, 
as we shall see later. In this case, the specific choice of the neutrino potential profile 
A(r) = y/2GfNe- turns out to have only a minor impact on the evolution of self- 
interaction effects, even if A(r) ^ f^{i^)- For the sake of definiteness, we single out one 
of the time-dependent profiles studied in [7] (the one at post-bounce time t = 5 s), and 
steepen it in a few km range above the neutrino sphere as suggested in [17, 50].§ The 
resulting profile A(r) adopted in this work is shown in Figure 2. 

§ The local steepening does not change any feature of the collective neutrino flavor transformations, 
but it turns out to help the start-up of our numerical evolution routines in the first few km. 
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200 



Figure 2. Radial profiles of the neutrino self- interaction parameter //(r) = 
\/2Gf (N + N) and of the matter-interaction parameter A(r) = \/2Gf N^- adopted 
in this work, in the range r G [10, 200] km. Shown are also the approximative ranges 
where self-interaction effects are expected to produce mainly synchronization, bipolar 
oscillations and spectral split (see the text for details). 



3. Single-Angle Approximation: Analytical aspects 

An ensemble of relativistic neutrinos and antineutrinos coming in F flavors can be 
described by a set of dimensionless F x F density matrices pp and pp, one for each 
momentum mode. The most general Liouville evolution equations for pp have been 
worked out in [51]. In this Section we specialize and discuss such equations in Bloch form 
for F = 2 and single-angle approximation. Although we mostly rely on the formalism 
introduced in [43] and on the results presented in [17,42-45], we think it useful — for the 
sake of clarity and self-consistency — to give a systematic overview of the asymmetric 
{ue 7^ Ue) and continuous-energy case, also because our choice of thermal spectra allows 
some useful analytical estimates in terms of temperature parameters (3 = 1/T. 
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3.1. Block Vector Notation 

We consider two-family mixing between Ue and ly^ {Ve ^-nd Vx), driven by mass- mixing 
parameters (Am^, sin^ ^13). In this case, one can switch from normal to inverted mass 
hierarchy in two ways: (1) mixing is kept unaltered while -|-Am^ — > — Am^; or (2) Am^ 
is kept positive, while the mixing angle is swapped, sin0i3 cos ^13 (which implies 
-|- cos 2^13 — > — cos 2^13, with unaltered sin 2^13). Hereafter we adopt the latter choice, 
as advocated in [43]. In numerical calculations, we assume reference values 

Am^ = 2 X 10-2 eV^ , sin^ = 10"^ . (18) 

In the single-angle approximation, the only kinematical parameter is E — |p|. 
For each energy E and radius r, the neutrino density matrix p in flavor basis can be 

written in terms of polarization (Bloch) vector P = (Pj,, P^,, Pj)^ = P^jX + Pj^y + P^z, 
being x, y and z three orthogonal unit vectors, and of the vector of Pauli matrices 



Pee Pe.\l \ f -eK \n^^_ ^^^^ 
V Pxe Pxx J V ^e^^o; Wx\ J 2 

where 1 is the unit matrix, n — Tr(p) = Ue + Ux represents the occupation number of 
neutrinos (per unit volume) at energy E, and the densities are defined through Eq. (11) 
in our model. Analogous definitions in terms of P hold for antineutrinos. The initial 
conditions in flavor basis p* = diag(ne, Ux) and = diag(ne. Tlx) correspond to 

pi^pir,^!!:!^^ (20) 

F = i?z=^^^z. (21) 

n 

The final (/) survival probabilities are then given by 

Pee = P{K - l^i) = \{l + ^^ , (22) 




Pee = P(l7:^I7/) = - 1 + ^ . (23) 



It is useful to introduce the integral polarization vectors of neutrinos and 
antineutrinos. 



J = 



L=jdEn^, (24) 

L_JdEnF, (25) 
as well as their sum S and difference D, 



J = 

N + N 



S = J + J , (26) 
D = J - J . (27) 
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The initial conditions imply that 

N + N /?, + + 2/3, ' ^ ' 

T^gl^z^ z, (29) 

N + N /3e + /5, + 2/5, ^ ' 

S'=^^i±^i^z = ^±i^z, (30) 
N + N /5e + /5, + 2/3, ^ ^ 

D'^^^z^ z. (31) 

iV + iV (3^ + ^^ + 2/3, ^ ^ 

Another auxiliary (unit) vector is the "magnetic field," 

B = sin 2^13 cos 26*13 z (32) 

where the upper (lower) sign refers to normal (inverted) hierarchy. For small ^13, one 
can take B ~ +z, unless the relevant dynamics is very close to the z axis. 

3.2. Equations Of Motion And Removal of Matter Effects 

The Bloch equations of motion (EOM) for each polarization vector P read, in single- 
angle approximation, 

P = i+ouB + Xz + //D) X P , (33) 
F = {-ouB + \z + //D) X P , (34) 

where the three terms in brackets embed vacuum, matter, and self-interaction effects. In 
particular, the third term couples all modes P and P, and is responsible for collective 
effects. The equations conserve each |P| and thus unitarity. It is understood that 
P = P{E,r), A = A(r), and /i — iJ>{r), with r — t and — dr- If the continuous 
parameter E is discretized through a set of Ne points {Eh}h=i,...,NE: then a set of 
6 X Ne coupled, first-order differential equations in t is obtained. The equations are 
"stiff," namely, their solutions generally involve a fast-changing combination of multi- 
frequency oscillations, due to the "precession" of P and P around the three terms in 
brackets. It is thus amazing that, through appropriate approximations, the resulting 
dynamics turns out to be understandable in terms of simple phenomena, as discussed 
below. 

In a frame rotating with angular velocity Az [42], the time derivative acquires an 
extra operator — Azx, which cancels the matter term. Moreover, the z-component 

of P (and of any other vector) are unchanged in such co-rotating frame, and 
thus the survival probability Pee is also unchanged. Only the transverse {x, y) 
components of any vector get mixed, e.g., the "magnetic field" becomes B = 
(sin 2^^i3 cos(— At), sin 2^^i3 sin(— At), =1= cos 2^^i3)"^. Again, one can still take B ~ =i=z, 
unless the relevant dynamics does not occur too close to B or z (a case which will be 
discussed separately). Matter (A) effects thus "disappear" in such co- rotating frame, as 
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pointed out in [42], leaving the other terms in the EOM formally unchanged: 

P = (+a;B + //D) X P , (35) 

F = (-a;B + /iD) X P . (36) 
The corresponding EOM for integral quantities are: 

j = + BxW + /iDxJ, (37) 

j = - BxW + /xDxJ, (38) 

S = B X (W - W) + //D X S , (39) 

D = Bx(W + W), (40) 



W 



where 

— ^ / dEcunP , (41) 

W= — ^-^fdEcunP, (42) 
with initial conditions implying 

= dEuj(ne-nJz^ ^ z, (43) 



W = = / dEuj (rie -nJz^ _ ^ z . (44) 

N + nJ ^ ' 2 B. + B. + 23^ ^ ^ 



5.5. Conservation Laws and Spectral Split 

The equation for D imphes conservation of the scalar 

D • B = const = • B ~ ^D* • z = +— — =J- , (45) 

corresponding to the conservation of the (electron) lepton number, and implying pair 
conversions of the kind i/gl^e t^x^x [43] . 

In the special case of (nearly) constant neutrino density (/i ~ 0), another conserved 
scalar is the average energy per neutrino pair, given by 

£: = B-(W + W) + ^//D2 = V + r , (46) 

where the first (second) term acts as a sort of potential (kinetic) energy. 

When r is sufficiently large to make self- interaction effects vanish (/x <C a;), the 
kinetic term T also vanishes, and the "ground state" for E would corresponds to the 
minimization of the potential V, i.e., to W and W aligned as much as possible with 
— B (in any hierarchy), provided that lepton number is also conserved. 

Since the vectors W and W always start aligned to z, this implies that, in normal 
hierarchy (z ~ — B) they end up in the same position, trivially conserving lepton number 
(i.e., the system starts — and remains — close to the minimum of the potential energy). 
Conversely, in inverted hierarchy (z ~ +B) the vectors W and W start antialigned 
with — B (maximum of the potential), and for large r they should reverse their direction 
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in order to approach the potential minimum. Reversal cannot be complete, however, 
since it would maximally violate (invert) lepton number. Consistent minimization of the 
potential can be achieved by complete reversal of the smallest vector W, and by partial 
reversal of the largest vector W (note that the excess of neutrinos over antineutrinos 
leads to |W| > |W|). 

More precisely, only a fraction W> of W, the fraction above a certain critical (split) 
energy Ec, is reversed, while the complementary fraction, W< = W — W>, remains 
unaltered [44-49]. In such final state, the critical energy is fixed by lepton number 
conservation, i.e., by Dl — D( with 



the last two terms having the opposite overall sign in the initial state {N -\- N)D\. Then 
one gets an implicit equation for Ec, 



which, in our specific SN model, is solved for Ec — 7 MeV. When all collective effects 
are terminated (// <^ cu), one expects a nearly complete inversion of the polarization 
vectors for E > Ec, and thus a "stepwise swap" between the and i/x energy spectra. 

Of course, such reasoning is heuristic, and does not prove that the dynamics allows 
the system to reach the peculiar final state discussed above. We refer the reader to [44- 
48] for explicit solutions constructed in adiabatic approximation (slowly decreasing /i), 
which indeed lead to spectral split under rather broad assumptions. On the other hand, 
such adiabatic solutions average out the interesting transient phenomenon of bipolar 
oscillations [42,43], which we discuss below. 

3-4- Alignment Approximation and Gyroscopic Pendulum 

For /i|D| ^ uj, Eqs. (35)-(38) reduce to the same form V ^ yuD x V, and thus 
all polarization vectors V = P, P, J, and J have the same dynamics (in particular, 
they remain closely aligned to each other, and to the 2;-axis, as they are at the start). 
As fjL decreases, the ±a; terms in the EOM start to be non-negligible, and neutrino 
and antineutrino polarization vectors develop different precession histories. As far as 
Acu/ii remains small, where Acu is the typical energy spread, the individual P's stick 
to the global vector J (i.e., their components parallel to J typically dominate over the 
transverse ones), and analogously the P's stick to J, with J and J gradually separating 
from each other. || 

Within such "alignment approximation," also W (W) is nearly parallel to J (J), 



II This approximation may become ill-defined in the symmetric case rie = Tie (not our case), where the 
smallness of |D| makes the condition /xjDj ;» co critical and the (J, J) separation unclear. See, e.g., the 
remarks in Sec. VI A of [43]. 





(48) 



W ~ w J , 
W ~ WJ , 



(49) 
(50) 
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and Eqs. (37)- (38) become 

j = [(^^dif + t^ave)B + ^D] X J , (51) 

J = [{uJdif - a;ave)B + /iD] x J , (52) 
where we have defined 

t^ave ^ (w + w)/2 , (53) 

t^dif = (w - w)/2 . (54) 

Equations (51)-(52) imply conservation of the vector moduli J = |J| and J = |J|, 
as well as of 1^ = |W| and W — |W|. The frequencies w and w can be evaluated, e.g., 
in the initial state, providing 

= 1 . = IdEujne-nx) JdE^jn^-fix) 

= ^^(/9e + /3, + 2/5.) . (55) 
In our specific SN model, for Am^ = 2 x 10~^ eV^ it is 

u;ave ^ 0.9 km~^ . (56) 

By going in a co-rotating frame with frequency a;difB ~ ^cUdifZ, the terms a;difBx in 
Eqs. (51)-(52) are rotated away, with no other formal change in the EOM for J and J. 

The a;difBx terms disappear also from the EOM of S and D. By defining 
Q = S — (a;ave/A*)B [43], and assuming /i ~ (adiabatic variations of //), Eqs. (39)-(40) 
can be written as 

Q = //D X Q , (57) 
D = a;aveB x Q , (58) 

showing that the (i^, I^) ensamble is characterized by a single, collective frequency cuave, 
despite the existence of a continuous energy spectra. Notice that Q — |Q| is conserved, 
as well as D • B and D • Q [43]. 

It has been realized that Eqs. (57)- (58) describe a gyroscopic pendulum in fiavor 
space [43,44]. In particular, by making the identifications 

Q/Q = r (unit length vector) , (59) 

D = L (total angular momentum) , (60) 

fj,~^ = m (mass) , (61) 

D • Q/g = a (spin) , (62) 

a;ave A* Q B = - g (gravity field) , (63) 

one can write Eqs. (57)- (58) in the form (after right-multiplying Eq. (57) by rx) 

L = mv XV -\- ar , (64) 
L = mr X g , (65) 
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which are the equations of motion of a spherical pendulum of unit length (|r| = 1), 
subject to a constant gravity field g, and characterized by a point-like bob of mass m 
which spins around the pendulum axis r with constant (inner) angular momentum cr. 
The most general evolution of this system is a combination of two periodic motions of 
the bob, one around the vertical g axis (precession) and the other along it (nutation) 
[43,44]. An exphcit solution of the EOM can be constructed in terms of quadratures 
involving the three integrals of motion [52, 53], namely, the spin a — L ■ r, the vertical 
component of the angular momentum L • g/|g|, and the energy £ which now equals 

£: = -mg-r+ -1-2+— , (66) 



where we have kept the spin term in the (bracketed) kinetic energy. The explicit 
solution, however, involves elliptic integrals and is not particularly transparent. Here it 
is sufficient to recall the following global features of the pendulum motion. 

In normal hierarchy, the pendulum starts close to the stable, downward position 
(the misalignment being of 0(^^13)), and remains close to it as slowly decreases (i.e., 
m slowly increases). Conversely, in inverted hierarchy, the pendulum starts close to 
the "unstable," upward position. When is large and thus m is small, however, the 
bob spin a dominates ("fast rotator"), and the pendulum remains precessing in the 
upward position to conserve angular momentum ("sleeping top" ) [44, 52]. This situation 
(also dubbed as "synchronization" [32, 43] in the SN neutrino context) is stable if the 
dominant (spin) kinetic term is larger than the maximum excursion of the potential 
energy [43,44,52], namely, 

> 2m g . 67 

2m ^ ' 

For large ixjoj^^^ (and thus vertically aligned polarization vectors), this condition 
translates into /i > /igup, where 

4awc^_, (7Ve + A^e)^-(iV. + A^.)^ 

In our reference SN model, it is 

//sup ~ 75 a;ave - 67 km~^ . (69) 

Conversely, when // < //gup, any initial misalignment with the "vertical axis" set 
by g (i.e., any ^13 ^ 0, no matter how small) triggers the first fall of the pendulum 
and its subsequent nutations (also dubbed as "bipolar oscillations" in the SN neutrino 
context). If n were constant, bipolar oscillations would be exactly periodic. However, 
the decrease of yU implies an increase of the pendulum inertia (m); the pendulum never 
swings back to exactly the same uppermost position, which instead steadily decreases, 
together with the vertical amplitude of nutations (roughly oc jd}^"^ [43,44]). Bipolar 
oscillations are then expected to vanish when self-interaction and vacuum effects are 
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comparable, and the "alignment approximation" breaks down. More precisely, one may 
expect this condition to occur when the and u terms in the EOM of J and J are of 
the same size (/iD • B ~ a;ave), implying fi ~ /Xinf with 

_ CUave _ N + N _ /3e+^e + 2/3^ 

/^inf - - '^ave^^^— ^ - '^ave _ (70) 

In our reference SN model, it is 

/Xinf ~ 7.5 u;ave ^ 6.7 km"^ . (71) 

The condition fj, ~ //inf roughly marks the "end" of bipolar oscillations and of 
collective effects, but not yet of all self-interaction effects. In particular, spectral split 
effects continue to build up for //<//inf, and eventually freeze out for /i <^ (^ave- The 
reason is that the neutrino spectral split requires a separation of the vector W into 
parts W< and W> oppositely evolving in the z-component. As far as the alignment 
approximation holds (and thus bipolar oscillations occur), this split cannot fully develop, 
and should therefore be complete somewhat beyond the bipolar range. Of course, there 
is no sharp boundary between the two processes: in the range where /i ~ /^inf , one should 
observe a smooth vanishing of bipolar oscillations, and a gradual build-up of the spectral 
split, through the polarization reversal of neutrinos with E > Ec, with an associated 
non-conservation of W (and of J). Summarizing, we expect the following sequence of 
dominant phenomena, where the radial ranges refer to our reference SN model: 

A* ^ A*sup ■ synchronized oscillations (r < 55 km) , (72) 

A*inf ^ A* ^ A*sup : bipolar oscillations (55 km < r < 100 km) , (73) 
fJ' ^ A^inf : spectral split (r > 100 km) . (74) 

Such ranges are explictly shown in Fig. 2. For numerical purposes, we shall stop our 
investigations to 200 km in this paper. We are not concerned here with subsequent 
(ordinary or "stochastic") MSW transitions which may occur later at r ~ O(IO^) km 
when A ~ a;. 



3.5. The Revenge of Matter Effects 

Matter A effects can alter the previous description in two ways: (1) by anticipating the 
MSW condition A ~ a;ave within the first 0(200) km in shallow matter profiles (not our 
case), thus interfering with collective effects in the same range [17, 46]; (2) by altering 
the dynamics when the polarization vectors are very close to z, which occurs just at the 
transition between synchronized and bipolar oscillations [43]. When the latter transition 
occurs, matter (A) and small mixing (^13) effects cannot be completely co-rotated away, 
the B ~ =l=z approximation fails, and the transverse components of B (oscillating with 
amplitude sin 26*13 and frequency A) must be taken into account. In the nontrivial case 
of inverted hierarchy, it turns out that their general effect is to further stabihze the 
"upward" pendulum position [43], elongating the period of (at least) the first bipolar 
swing. Explicit analytic estimates have been presented in [43] for the symmetric case 
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Figure 3. Single-angle simulation in inverted hierarchy: Modulus and z-component 
of the global polarization vector of neutrinos (J) and antincutrinos (J), as a function 
of radius. The difference = Jz — Jz (not shown) remains strictly constant in r. 



rig = He, by solving the full EOM with time- varying B and for small deviations around 
the vertical directions. If we take these estimates as a reasonable proxy also for our 
asymmetric case {ug ^ n^, we expect that the first bipolar oscillation (nutation) should 
be a factor of ~ 4 longer than the next ones (which are much less affected by these 
subtle matter effect, the uppermost pendulum position becoming increasingly tilted 
for decreasing /i). The onset of bipolar oscillations is then expected to be somewhat 
delayed beyond r = 55 km in our reference SN model, by a time span equivalent to "a 
few nutations." 

4. Single-Angle Approximation: Numerical results 

The previous analytical expectations for the single-angle case are nicely confirmed by our 
simulations. We numerically evolve Eqs. (33)- (34) in the range r G [10, 200] km within 
our reference SN model, considering only in the nontrivial case of inverted hierarchy (we 
have anyway checked that no significant effect occurs in normal hierarchy). Technical 
details are discussed in the Appendix: here we focus only on the results and their 
interpretation. 

Figure 3 shows the radial evolution of the modulus J and ^-component of the 
global neutrino polarization vector J (and analogously for the antineutrino vector J). 
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Their radial profile can be interpreted as follows. Up to ~ 68 km, it is J = and 
J = J z'- all polarization vectors are "glued" (synchronized) along the vertical axis, 
and the gyroscopic flavor pendulum just spins in the upward position without falling. 
At r ~ 68 km, the pendulum falls for the first time and nutations appear, marking 
the transition from synchronized to bipolar regime. The transition is retarded by a few 
nutation periods (from the expected ~ 55 km to ~ 68 km) by the matter effects discussed 
in the previous Section. The nutation amplitude gradually decrease (as ~ A*^^^), and 
bipolar oscillations eventually vanish for r > 100 km, as expected. 

At the same time, the spectral split effect builds up. Antineutrinos tend 
to completely reverse the polarization vector (J thus minimizing their 

"potential energy" (after which nothing relevant happens to them), so that — — J 
asymptotically.^ Neutrinos also try to invert their global polarization vector (as much 
as it is allowed by lepton number conservation) as soon as the alignment approximation 
breaks down (Ai</imf) and non-conservation of J is allowed. Indeed, for r>100 km, 
J decreases. Eventually the situation J ~ is reached (slightly beyond the r range 
in Fig. 3), when the spectral splitting is frozen, corresponding to a final state with J 
aligned with -\-z and J aligned with — z. In all the above processes, lepton number is 
strictly conserved, leading to the constancy of D^ — Jz — J z at any r. Prom the point of 
view of observable oscillation probabilities (related to the z-component of polarization 
vectors), the situation is basically frozen well before r ~ 200 km. We can conclude 
that the behavior of J, J, and J ^ is well understood, with good agreement between 
analytical expectations and numerical simulations. The agreement can also be extended 
to more detailed features of the energy spectrum, as discussed next. 

Figure 4 shows the behavior of the individual polarization components Pz (upper 
panel) and Pz (lower panel) as a function of r, for five representative values of energy (in 
MeV): Ex ~ 5.2, ~ 12.4, E'i ~ 19.1, E^ ~ 23.8, and E^ ~ 31.3 (which are a subset of 
the grid sampling energies, hence their "non-rounded" values)."*" In Fig. 4, the onset of 
bipolar oscillations and their nutation periods are clearly equal for both v and V at any 
energy, confirming the appearance of a self-induced collective behavior, characterized by 
a single frequency parameter cjavc for all (anti)neutrinos, despite the spread of vacuum 
oscillation frequencies uj. However, the fate of each Pz or Pz does depend on energy. For 
neutrinos (upper panel in Fig. 4), one can clearly see the phenomenon of spectral split 
around -Ec — 7 MeV: the curve at Ex < E^ ends up at the same initial value, while the 
curves at £^2, £'4, -Bs > Ec show the expected inversion Pz — > —Pz- Only the curve at 
energy £'3 does not change much {Pz ~ at any r) , being close to the energy where 
the Vf. and fiuxes are equal (see Fig. 1), and fiavor transformations are inoperative. 
For antineutrinos (lower panel in Fig. 4), all curves show complete polarization reversal 
as expected [P z ~^ —Pz), including the "trivial" case ~ at £4 ^ E^q, where the 
Ve and fiuxes are equal. We conclude that the numerical simulations confirm the 

^ There is actually a slight loss of J during bipolar oscillations in Fig. 3 (as also numerically observed 
in [49]) which makes the final |J/| ~ 0.09 slightly smaller than the initial J* ~ 0.1. 
"•" We do not show the moduli, which are always strictly conserved. 
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Figure 4. Single-angle simulation in inverted hierarchy: z-component of the 
polarization vector of neutrinos (P^ , upper panel) and antineutrinos (P^ , lower panel) 
as a function of radius, for five representative values of the energy. 



end of bipolar oscillations and the appearance of the energy split phenomenon (around 
r ~ 100 km) with the expected global features. There is only a minor "unexpected" 
effect (lack of Pz reversal for E <4: MeV, not shown in Fig. 4), as commented below. 

At r = 200 km, self-induced flavor transformations have basically ended. The 
transition from initial (r = 10 km) to final (r = 200 km) fluxes implies, in our two- 
family scenario (Pge = Pxx = 1 — Pex, see also the remarks at the beginning of Sec. 2), 

^ mp,^^E),m^-p.AE)], (75) 



(B,) (B,) ' {E,} 

and similarly for antineutrinos 
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Figure 5. Single-angle simulation in inverted hierarchy: Final fluxes (at r = 200 km, 
in arbitrary units) for different neutrino species as a function of energy. Initial fluxes 
are shown as dotted lines to guide the eye. 



Figure 5 shows the final neutrino and antineutrino fiuxes calculated in this way. 
The left panel (neutrinos) clearly shows the spectral split effect, and the corresponding 
sudden swap of z/g and fluxes above Ec — 7 MeV. In the right panel of Fig. 5, 
the flnal antineutrino spectra are basically completely swapped with respect to the 
initial ones (compare with Fig. 1), except at very low energies, where there appears 
an "antineutrino" spectral split. We relate this phenomenon to the loss of J and 
of \Jz\ observed and commented in Fig. 3: The small deflcit < X, can indeed 
obtained, analogously to the neutrino case, through the lack of Pz reversal for low-energy 
antineutrinos {E < E^ with Ec<4: MeV). A better understanding of this minor feature 
and of the \Jz\ loss would be desirable; however, we anticipate that the antineutrino 
spectral split is largely smeared out in multi-angle simulation, contrary to the neutrino 
spectral split which appears to be a robust, observable feature. 



5. Multi-angle Simulations: Notation and Numerical results 



In the multi-angle case (applied to the neutrino bulb model), any single polarization 
vector depends not only on the energy E and on the total propagation distance along a 
neutrino trajectory t, but also on the (incident) angle between the neutrino trajectory 
and the polar axis. The vectors Pi){E,t) and P^{E,t) obey then the following EOM 
[17], 

P^= +uB+Xz+2nV2GFjdc^>dE{l-c^c^>){jP./,>-JPi,>) xP^,(77) 
Ftf = -uB+Xz+27rV2GF fdc^>dE{l-c^c^^){jP^> -jP^>) xP^,(78) 
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where = COST? with ■& e [0,'i?max] and c^/ = cos'd' G [cos'i^max, 1],* while j and j are 
the total neutrino and antineutrino number densities 

j = je + jx , (79) 

J-Je + Jx, (80) 

as defined in Sec. 2.2. 

It is useful to characterize all the properties of the neutrino beam by using the 
radius r and the emission angle ■j?o at the neutrinosphere [17], by means of 

r sm'& — R,y siw&o , (81) 



'r^-Rlsin^^o-RuCos^o , (82) 

the range of cos 'do being constantly [0, 1] at any r. Along a generic trajectory at angle i), 
dt = dr/c^, and thus dt = c^dr. Since any couple (t,"*?) is in one-to-one correspondence 
with (r, t^o)) the polarization vectors can be relabeled as P^g(i?,r). 

For convenience, one may define effective densities n as in the single-angle case, 
n(r, E) = 27rD{r)j{E), so that the initial conditions ai r — Ri, become as usual 

= Ue - Jx)/j = {ue -n^)/n , (83) 

P'tfo = Ge - D/J = - n^)/n . (84) 

The single-angle limit is recovered by fixing d = Q = {}q and by assuming that 
all polarization vectors behave as the ones at -i^o = 0, in which case the integral 
/ dc^i{l — c^i) = D{r) is factorized out (this is actually the way D{r) is originally 
defined [17]). ^ 

Equations (77)-(78) reduce then to Eqs. (33)-(34) after the integrated densities N 
and N are introduced as in Sec. 2.2. 

Angle-averaged polarization vectors can be defined as 

P = , (85) 

(and similarly for P), where the "extra" cosine factor accounts for projection in radial 
direction [49]. Global polarization vectors J, J and D = J — J can then be defined in the 
same way as in Sec. 3.2. Although their EOM are not as simple as in the single-angle 
case, it turns out that D ■ B ~ is still a conserved scalar [43,49] (proof omitted). 

The angular dependence of the neutrino-neutrino interaction strength is generally 
expected to introduce some "dephasing" or kinematical decoherence between different 
neutrino trajectories, and to smear out the "fine structures" observed in single-angle 
simulations. While decoherence effects are dominant in the symmetric case {rie = fie) 
[34], they seem to be only subdominant in asymmetric cases {rig ^ He) [17,49]. In 
the latter case their quantitative description lacks, at the moment, of an analytical 

* The angles and '&' must lie in the cone subtending the neutrinosphere. The factor 2-k = J dip comes 
from cylindrical symmetry within this cone. 

tt An alternative single- angle case has been recently studied in [49], by selecting the emission angle 
'&0 = 7r/4 instead of do = = 
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Figure 6. Multi-angle simulation in inverted hierarchy: Modulus and z-component of 
the global polarization vector of neutrinos (J) and antineutrinos (J), as a function of 
radius. The difference = Jz — J z (not shown) remains strictly constant in r. 

understanding, and relies mainly on numerical simulations [49]. However, one may 
at least expect that the appearance of the neutrino spectral split phenomenon is 
not spoiled in multi-angle cases, being based on the broad-brush picture of potential 
energy minimization (i.e., final (anti) alignment with B) constrained by lepton number 
conservation (i.e., constant D^)- This expectation is confirmed by the numerical 
simulations discussed below, which refer only to the nontrivial case of inverted hierarchy. 

Figure 6 is the multi- angle analogue of Fig. 3. By comparing the two figures, 
it appears that bipolar oscillations of J and J are largely suppressed in the multi- 
angle case, only the first nutation being clearly visible. Moreover, such nutation starts 
somewhat later (r ~ 76 km) as compared with the single-angle case (r ~ 68 km). These 
features can be understood in terms of the different self-interaction effects experienced 
along different trajectories. In multi-angle simulations, neutrino-neutrino angles can be 
larger than the (single-angle) average one, leading to somewhat stronger self-interaction 
effects, which keep the system in synchronized mode for a slightly longer time, and thus 
delay the first nutation. Along different trajectories, the subsequent bipolar oscillations 
have also somewhat different amplitudes and phases, which tend to cancel out in the 
global polarization vectors. For a similar reason (relatively stronger self-interaction 
effects, as compared to single- angle), in Fig. 6 there is a slightly more pronounced loss 
(depolarization) of J at the start of the first bipolar oscillation, and at the same time a 
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Figure 7. Multi-angle simulation in inverted hierarchy: z-component of the 
polarization vector of neutrinos (P^ , upper panel) and antineutrinos (P^ , lower panel) 
as a function of radius, for five representative values of the energy. 

prolonged "coherence" of J (slower decrease of J), with respect to Fig. 3. However, just 
as in the single-angle case, it turns out that gets finally reversed, while the difference 
Dz = Jz — J z is exactly conserved. The reversal becomes more evident by looking at 
specific (anti)neutrino energies. 

Figure 7 shows the behavior of the z-component of angle-averaged polarization 
vectors [Eq. (85)] for neutrinos (P^, upper panel) and antineutrinos (P^, lower panel) 
in our multi-angle simulation, at fixed energies. The behavior is qualitatively similar 
to a "smeared version" of the single-angle curves in Fig. 4, with polarization vectors 
reversing (or not) their z-components as expected. The small depth (or absence) of 
nutations makes it more evident that the polarization reversal (i.e., the spectral split) 
starts to dominate over the bipolar mode around the expected radius r ~ 100 km. 
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Figure 8. Multi-angle simulation in inverted hierarchy: Final fluxes (at r — 200 km, 
in arbitrary units) for different neutrino species as a function of energy. Initial fluxes 
are shown as dotted lines to guide the eye. 

Figure 8 shows the final (r = 200 km) z/ and V fiuxes as a function of energy. The 
neutrino spectral swap at E > E^. — 7 MeV is rather evident in the left panel, although 
it is less sharp with respect to the single-angle case in Fig. 5. In the right panel of Fig. 8, 
the minor feature associated to the "antineutrino spectral split" is largely smeared out 
(see the same panel in Fig. 5), and survives as a small excess of at low energy. 

The spectra in Figure 8 are largely independent from the specific mixing value 
chosen for the simulations (sin^ 6*13 = 10~^), as far as 6*13 > (as we have also 
checked numerically). Variations of sin^ ^13 only lead to logarithmic variations in 
the (unobservable) synchronized-bipolar transition radius, and in the depth of bipolar 
oscillations [43,44], which are anyway smeared out in multi-angle simulations, as we 
have just seen. Therefore, the spectra in Figure 8 may be taken as rather general 
"initial conditions" for possible later (ordinary or stochastic) matter effects, occurring 
when ~ A(r) at r ^ 200 km. These later, ordinary matter effects are instead strongly 
dependent on 613, and vanish for, say, sin^ 6*13 < 10~^ (see, e.g., [7]). If 6*13 is indeed that 
small (but nonzero), neutrino self- interaction effects could be the only source of fiavor 
transformations in (anti) neutrino spectra. 

In conclusion, for < sin^ 6*13 < 10~^, the observable spectra at the SN exit 
would be similar to those in Fig. 1 for the normal hierarchy case (no significant fiavor 
transformations of any kind), and to those in Fig. 8 for the inverted hierarchy case (large 
self-interaction effects). For sin^ ^13 > 10~^, the same spectra should be taken as "initial 
conditions" for the calculation of subsequent MSW effects. Once more, we remark that 
the decouphng of self-interaction and MSW effects is a characteristic of our adopted 
SN model, inspired by shock- wave simulations [7]. The phenomenology becomes more 
complicated in alternative models with shallow matter profiles, when both effects can 
occur in the same region, as in the simulations performed in [17,47]. 
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6. Conclusions and prospects 

Neutrino-neutrino interactions in the high-density region near a supernova core have 
been recently recognized to produce surprising collective effects. Motivated by these 
developments, we have investigated neutrino flavor transformations in a SN model, 
where the main self-interaction effects (synchronization, bipolar oscillations, and 
spectral spht) develop well before possible MSW effects. 

The neutrino-neutrino interaction strength depends on the intersection angle of 
their trajectories. Averaging this (variable) angle along a single, representative radial 
trajectory leads to the so-called single-angle approximation, which allows both elegant 
analytical insights [43,44] and easier numerical calculations. However, removal of this 
approximation is needed (through multi-angle simulations) in order to validate the 
analytical insights, and to show that the main effects are not spoiled by kinematical 
decoherence. Moreover, many self-interaction effects have actually been first seen 
numerically and then interpreted analytically a posteriori. 

We have thus performed numerical simulations in both single- and multi-angle cases, 
using continuous energy spectra with significant I'-T' and Ue-i^x asymmetry. The single- 
angle results can be understood analytically to a large extent, and their main observable 
effect — in the nontrivial case of inverted hierarchy — is the swap of energy spectra above a 
critical energy dictated by lepton number conservation [45]. In multi-angle simulations, 
we find that the "fine structure" of self- interaction effects (e.g., bipolar oscillations) is 
smeared out, but the spectral swap remains a robust, observable feature. In this sense, 
trajectory averaging does not play a crucial role. This is the main result of our work. 

The swapping of the and Vfj^ (as well as of the and i^^) fluxes could have an 
impact on r- process nucleosynthesis [35-37, 50] , on the energy transfer to the stalling 
shock wave [54], and on the possibility to observe shock- wave propagation effects in 
neutrinos. It is also worth studying possible self-interaction effects in the phenomenology 
of SN 1987A neutrino events [55-57] and of the diffuse supernova neutrino background 
spectrum [58,59]. Further analytical and numerical developments may require to solve 
the neutrino evolution equations in the complete Su flavor scenario, where new effects 
associated to the "solar" Sni^ can occur; a recent example has been worked out in [60]. 
Perturbations of the (cylindrically symmetric) bulb model for neutrino emission might 
also be considered in more advanced simulations. Finally, there is a continuing interest 
in more formal aspects of the mean-field approach to neutrino self-interaction effects 
[61], which has been implicitly assumed in most of the related literature (including this 
work); see [62] for a recent discussion of its validity and inherent approximations. 

In conclusion, twenty years after the SN1987A, the understanding of SN u flavor 
transformations is still in progress, and surprising self-interaction effects are emerging as 
possible dominant phenomena. These effects are changing the current paradigm of SN 
neutrino physics, and demand further analytical and numerical investigations, as well 
as new experimental inputs and guidance — should Nature be so kind to make a galactic 
supernova explode. 
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Figure 9. Absolute difference between numerical evaluations of Jz{r) in multi-angle 
simulations, using various (energy) x (angle) grid sizes: 32 x 80 (baseline), 32 x 60 
(dashed), and 24 x 80 (dotted). 
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Appendix 

We discuss here a few technical aspects of multi-angle numerical simulations, which 
are much more challenging than single-angle ones. Equations (77,78), after the proper 
transformation ^ ^ 'do, provide a set of 6 x Ne x N^^ ordinary differential equations 
(ODE) in r (in the real domain), where Ne and N^^ are the number of points sampling 
the (anti) neutrino energy E and emission angle cosine cosi?o, respectively. 

This ODE set is stiff, and needs appropriate routines to be solved numerically. 
After a careful comparison of public routines, we have adopted the GAMD software 
[63] , implemented on a A^^; x = 32 x 80 grid. Denser sampling in cos i^o is required, 
since the polarization vectors generally vary more rapidly in than in E. An exception 
would be provided by MSW effects interfering with self- interaction ones [17] (not our 
case), which requires much denser sampling, especially in E, in order to track the MSW 
resonant behavior. In our scenario with no MSW interference, we obtain satisfactory 
numerical convergence with Ne x N^^ = 32 x 80. Figure 9 compares the last three 
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steps in a trial sequence with increasingly denser grids A^^ x N^^ , showing that the final 
absolute error on the reference quantity can be safely estimated to be < 10^^. 

The grid points are not chosen to be equally spaced, but are instead fixed by 
imposing weighted Gaussian quadrature of the double integrals in the right-hand-side 
of Eqs. (77,78). In other words, the Ne x N^^ grid points not only sample the energy 
and angular evolution of the polarization vectors in the ODE set, but are also used to 
perform the inner Gaussian integrations at each evolutionary step. This "trick" saves a 
lot of computer time. Nevertheless, a typical simulation over a 32 x 80 grid takes ~ 400 
hours on our local computer facility (a cluster of four, 64-bit and 4Gb RAM processors 
at 2.4 GHz, with Fortran 90 codes running on a Linux platform). We plan to use more 
powerful (remote) facilities in future works on the subject. 
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